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2.  Objectives 


The  following  were  the  principal  objectives  of  this  project. 

1.  Modeling  of  multistructures  via  hp  interface  mixed  methods  for  do¬ 
main  decomposition  and  concatenation.  This  work  complemented  an 
ongoing  implementation  in  the  commercial  program  MSC-NASTRAN. 

2.  Development  of  hp  methods  that  are  robust  in  the  presence  of  locking 
and  boundary  layers.  In  particular,  development  of  mixed  hp  methods. 
Applications  to  plates,  shells,  non-linear  elasticity  and  non-Newtonian 
fluid  flow. 

3.  Analysis  of  hp  space  enrichment  methods  and  applications.  The  driving 
engineering  application  was  one  of  non-linear  contact  problems. 

4.  Investigation  of  the  flnite  element  approximation  of  spectra  of  non¬ 
compact  operators.  Application  in  non-linear  buckling  analysis.  (Work 
in  progress.) 


1 


3.  Summary  of  Accomplishments/Status  of  Effort 

1.  The  search  for  a  non-conforming  hp  method  that  can  be  used  for  domain 
decomposition  and  concatenation  was  completed  by  establishing  that  the 
hp  mortar  finite  element  method  is  essentially  stable  and  optimal  in  terms  of 
both  h  and  p  for  a  wide  variety  of  meshes.  These  results  were  communicated 
to  MSC-NASTRAN  developers  who  incorporated  a  related  hp  method  in 
their  code.  Remaining  work  (in  progress)  on  this  topic  consists  of  analyzing 
some  new,  simpler  versions  of  the  method  that  we  have  formulated,  including 
one  that  extends  to  three  dimensions. 

2.  Research  was  completed  on  mixed  hp  finite  element  methods  that  were 
developed  over  straight-sided  and  curved  quadrilaterals  and  triangles.  This 
research  will  form  the  building  blocks  for  future  applications  to  various  prob¬ 
lems  (e.g.  viscoelasticity).  Already,  these  methods  were  applied  to  a  non¬ 
linear  large  displacement  problem  and  shown  to  be  more  robust  than  the 
standard  formulation.  In  a  related  project,  stable  hp  mixed  methods  were 
developed  for  non-Newtonian  fluid  flow  such  as  the  upper  convected  Maxwell 
model. 

3.  The  investigation  of  hp  finite  element  methods  for  modeling  boundary 
layers  in  thin  domains  was  completed.  In  particular,  rules  for  mesh-degree 
selection  for  the  robust  analysis  of  plates  and  shells  with  corners  were  devel¬ 
oped,  to  take  into  account  both  the  effect  of  corner  singularities  and  boundary 
layers.  Our  prescriptions  for  handling  boundary  layers  have  been  incorpo¬ 
rated  into  e.g.  the  STRESS  CHECK  user’s  manual. 

4.  Motivated  by  an  engineering  solution  method  proposed  by  Barna  Szabo 
for  non-linear  contact  problems,  an  hp  space  enrichment  method  and  its 
applications  were  investigated.  Theorems  that  mathematically  validated  the 
method  were  established,  setting  the  framework  for  future  work  in  this  area. 

5.  Research  on  low-order  elements  for  plate  and  shell  problems  was  continued 
with  the  completion  of  two  papers  that  analyze  the  robustness  of  mixed  FE 
methods  in  terms  of  both  locking  and  boundary  layers. 

6.  A  method  for  the  linearized  analysis  of  non-linear  buckling  problems  pro¬ 
posed  by  Barna  Szabo  (for  STRESS  CHECK)  was  investigated  theoretically. 
Preliminary  results  on  mathematical  validity  and  finite  element  approxima¬ 
tion  were  obtained.  This  investigation  will  form  a  significant  topic  in  the 
research  to  be  carried  out  under  a  subsequent  AFOSR  grant. 
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4.  Accomplishments 


1.  Modeling  of  multistructures 

We  mathematically  and  computationally  investigated  non-conforming  hp 
methods  to  model  problems  over  a  domain  consisting  of  several  subdomains. 
The  motivations  were  to  (1)  distribute  the  discretization  of  a  complicated 
domain  among  several  users  by  breaking  it  into  subdomains  which  can  be 
recomposed  for  the  final  analysis,  (2)  allow  selective  refinement  in  regions 
where  it  is  needed,  by  extracting  the  region,  re-meshing  it  and  putting  it 
back,  (3)  discretize  different,  inter-connected  PDEs  over  adjoining  domains. 
Such  methods  also  have  applications  in  domain  decomposition  methods,  e.g. 
the  overlapping  Schwarz  algorithm. 

One  such  method,  involving  three  separate  fields,  was  developed  and 
tested  by  researchers  at  NASA  Langley  for  low-order  elements  (p=l  or 
2).  A  version  of  this  was  under  implementation  by  the  MacNeal-Schwendler 
Corporation  for  their  commercial  program  MSC-NASTRAN.  Since  MSC- 
NASTRAN  uses  high-order  elements  as  well  as  low-order  elements,  their 
requirement  was  that  the  method  used  be  such  that  it  is  stable  and  opti¬ 
mal  in  terms  of  all  h,p  and  hp  methods  that  the  user  may  choose  (including 
meshes  that  are  highly  graded).  The  question  then  became  one  of  developing 
a  non-conforming  method  that  satisfied  all  these  criteria. 

We  fully  solved  this  problem  (for  the  two-dimensional  case)  under  this 
grant,  obtaining  sharp  theoretical  and  computational  results.  These  were 
communicated  to  MacNeal-Schwendler  and  were  used  by  them  to  select  and 
fine-tune  the  final  form  of  the  method  they  implemented. 

More  specifically,  we  analyzed  several  hp  non-conforming  methods,  and 
were  able  to  show  in  [1],  that  the  so-called  “mortar  finite  element  method” 
is  an  ideal  choice.  This  is  one  of  several  methods  that  is  based  on  the  use 
of  a  Lagrange  multiplier.  If  ui  and  U2  are  the  test  or  trial  functions  used 
in  different  sub-domains  fli  and  Q2,  then  instead  of  the  usual  continuity 
condition 

ui  —  U2  on  r  =  rii  n 
the  mortar  FEM  uses  instead 

J  {ui  —  U2)X  ds  =  0  VA  G  5,  (1) 
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where  S'  is  a  finite  element  space  of  Lagrange  multipliers.  Note  that  the 
meshes  on  F  from  the  two  domains  and  ^2  will  not  in  general  coincide.  S 
is  a  space  of  piecewise  polynomials  on  the  mesh  on  F  induced  by  the  partition 
of  one  of  the  sides,  say  Qi.  If  piecewise  polynomials  of  degree  p  are  used  for 
ui  and  U2,  then  for  S  this  method  uses  piecewise  polynomials  of  degree  p 
again,  except  that  in  the  first  and  the  last  interval  of  the  mesh,  the  degree  is 
one  less,  i.e.  p  —  1  (see  [1]  for  details). 

In  the  past,  this  method  was  only  proposed  for  h  refinement,  and  that 
too  over  meshes  that  were  quasiuniform.  Through  theoretical  analysis  as  well 
as  computational  testing,  we  established  stability  and  optimality  results  for 
essentially  arbitrary  meshes,  and  with  respect  to  h,p  as  well  as  hp  refinement. 
In  particular,  exponential  convergence  was  established  for  this  method  when 
geometric  meshes  with  increasing  p  are  used. 

In  this  regard,  the  following  theorems  were  proved  for  a  model  elliptic 
problem  on  a  domain  Q  decomposed  into  subdomains  Q,i  [1]. 

Theorem  1.  Let  u  be  the  exact  solution,  satisfying  u  €  >  3/2(1  > 

7/4  if  p  varies),  and  Uh^p  its  hp  mortar  finite  element  approximation,  using 
piecewise  polynomials  of  degree  p  on  quasiuniform  meshes  {T^}  on  each  fij. 
Then 

where  p  =  min{/,p  +  1}  and  C  is  a  constant  independent  of  h,p  and  u. 

(Here,  ||  •  |kd  is  the  usual  component- wise  discrete  H^{Q)  norm  introduced 
by  the  non-conformity.) 

Theorem  2.  Let  ujv  denote  the  solution  of  the  mortar  hp  method  with  N 
degrees  of  freedom,  where  geometric  meshes  are  used  so  that  the  degree  p 
is  proportional  to  the  number  of  layers  n  around  each  corner  of  in  each 
sub-mesh.  Then  the  error  decreases  exponentially,  i.e.  there  exists  a  7  >  0 
such  that 

Ik  -  Uh,p\\i,d  < 


The  paper  [2]  contains  some  additional  results  that  indicate  that  the 
estimates  obtained  in  [1]  are  sharp,  and  also  that  give  computational  results 


for  a  different  test  problem.  We  found  that  characteristic  hp  convergence 
curves  are  obtained  (as  predicted  by  Theorem  2  above),  showing  very  little 
difference  from  the  conforming  method.  This  is  illustrated  in  Figure  1,  which 
shows  the  results  for  the  mortar  hp  method  over  an  L-shaped  domain  broken 
into  two  sub-domains.  The  mesh  is  geometrically  refined  over  each  domain, 
but  with  different  geometric  ratios,  so  that  the  meshes  do  not  match  on  the 
interface. 

These  results  will  be  further  augmented  in  [3].  There,  we  have  formulated 
a  much  simpler  version  of  the  mortar  method,  for  which  we  can  show  the  same 
numerical  performance,  and  establish  the  same  theoretical  estimates  as  the 
method  in  [1].  The  simplified  method  just  uses  degree  p  —  I  uniformly  over 
all  intervals  for  the  space  S.  (In  fact,  computational  results  indicate  that  the 
method  keeps  performing  well  even  when  degree  p  —  2  is  taken.) 

The  true  use  of  this  simplification  will  be  seen  in  three  dimensions.  This  is 
because  in  3-d,  the  mortar  finite  element  method  is  as  yet  not  well-formulated 
except  for  the  case  of  linear  elements.  Our  idea  of  taking  degree  k  —  1,  how¬ 
ever,  generalizes  easily  to  3-d.  The  convergence  properties  again  hinge  on  the 
stability  of  a  certain  projection  operator  11,  defined  essentially  by  equation 
(1)  (see  [1]).  Unlike  the  2-d  case  where  analytic  tools  were  used  to  character¬ 
ize  the  stability  of  H,  we  will  now  have  to  rely  on  some  computational  tests, 
which  we  are  in  the  process  of  developing.  These  results  will  once  again  be 
communicated  to  Dr.  John  Scheirmeier,  our  contact  at  MacNeal-Schwendler. 
This  3-d  work  will  be  the  final  section  in  the  thesis  [3],  which  has  an  expected 
completion  date  of  July,  1998. 

The  impact  of  our  results  will  be  mainly  via  incorporation  into  MSC- 
NASTRAN,  which  is  used  in  several  DOD  labs  and  aircraft  manufacturing 
companies.  There  are  several  3-d  components  (such  as  aircraft  parts)  for 
which  good  meshes  have  been  pre-computed.  This  technology  enables  users 
to  incorporate  them  in  the  final  mesh,  without  having  to  worry  about  in¬ 
terelement  continuity.  Moreover,  selective  refinement  of  the  mesh  or  the 
degree  can  be  implemented,  where  required. 

2.  Development  of  mixed  hp  finite  elements 

One  of  our  ongoing  goals  in  this  project  has  been  to  develop  and  analyze 
a  library  of  mixed  p  and  hp  finite  element  methods.  Such  methods  have  been 
very  useful  in  the  traditional  h  version  context,  where  they  have  been  used 
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in  problems  such  as  elasticity,  plates  and  shells,  where  they  overcome  the 
effects  of  locking.  They  have  also  proved  essential  in  other  applications  (e.g. 
Stokes  flow,  non-Newtonian  flow,  plasticity,  viscoelasticity  and  viscoplastic¬ 
ity),  where  the  variational  formulation  involves  a  saddle  point,  necessitating 
a  mixed  method.  It  is  anticipated  that  with  the  increasing  emphasis  on  p  and 
hp  methods,  and  the  application  of  these  methods  to  various  new  problems, 
a  library  of  validated  mixed  hp  elements  will  be  valuable  for  Air  Force  needs. 
More  specifically,  such  elements  are  essential  when  problems  involving  the 
constraint  of  incompressibility  are  to  be  solved,  such  as  e.g.  the  deformation 
of  rubber  tyres.  Future  applications  of  this  work  will  include  e.g.  the  fail¬ 
ure  of  electronic  components,  for  which  the  soldering  material  can  display 
viscoelastic  (incompressible)  behavior  at  high  temperatures. 

In  previous  work,  we  developed  and  analyzed  mixed  hp  methods  on  par¬ 
allelograms,  and  showed  that  with  these,  the  stresses  can  be  accurately  re¬ 
covered,  for  any  Poisson  ratio  (note  that  the  p  version  is  not  free  of  locking 
for  the  stresses).  Under  this  grant,  we  computationally  tested  these  ele¬ 
ments  and  also  extended  our  work  to  triangles  and  to  curved  elements.  This 
was  the  Ph.D.  dissertation  [4]  by  our  student,  Lt.  Col.  Lawrence  Chilton. 
Among  other  results,  we  (1)  developed  a  curved  quadrilateral  element  for 
which  we  can  theoretically  establish  the  same  stability  as  the  straight-sided 
parallelogram  element,  (2)  developed  mixed  hp  triangular  elements  which  are 
stable  in  h  and  demonstrate  an  inf-sup  constant  of  0{p~^)  in  the  range  of  p 
used  in  practice,  (3)  demonstrated  the  accuracy  of  these  elements  through 
computational  experiments. 

Let  us  describe  some  of  these  results  from  [4]  in  more  detail.  First,  for 
parallelograms,  we  considered  various  combinations,  the  most  well-known  of 
which  are  the  following.  (Since  these  mixed  elements  are  used  for  the  elas¬ 
ticity  equations  as  well  as  Stokes’  problem,  we  describe  them  for  the  latter, 
in  which  case  the  two  unknowns  are  the  velocity  and  the  pressure,  instead  of 
the  displacement  and  the  sum  of  the  normal  stresses.) 

(A)  [<3,1V<3p-2 

Here,  the  set  of  degrees  of  freedom  used  on  each  parallelogram  element  is 
Qp  for  each  component  of  the  velocity,  and  Qp-2  for  the  pressure.  This  el¬ 
ement  is  popular  especially  in  spectral  element  methods.  The  combination 
is  stable  in  terms  of  h,  while  the  stability  behaves  like  0{p~^^^)  in  terms  of 
the  degree.  .  In  terms  of  p,  the  element  gives  optimal  convergence  in  the  ve- 
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locity  (or  displacement)  and  no  worse  than  a  loss  of  in  the  pressure 

(or  stresses).  However,  it  does  not  give  the  optimal  approximation  of  the 
velocity  in  h,  since  the  pressure  degree  is  not  balanced  properly  according  to 
approximation  theory.  Hence,  it  is  not  recommended  for  hp  FEM. 

(B)  [QpIV.Pp-1 

This  combination  avoids  the  pitfalls  of  the  above  element,  by  being  optimal 
(and  stable)  in  h.  Recently,  it  was  shown  by  Bernardi  and  Maday  that  the 
stability  holds  in  p  as  well,  with  no  loss,  so  that  this  is  an  ideal  hp  element. 

(C)  [(3,nPp+2lVPp-i 

This  element  has  fewer  degrees  of  freedom  than  (B).  Although  this  reduction 
appears  to  reduce  the  theoretical  stability  estimate  to  0{p~^/'^),  our  compu¬ 
tational  results  in  [4,5]  indicate  that  it  behaves  as  well  as  (A). 

Let  us  now  describe  the  corresponding  methods  for  curved  elements, 
which  are  of  particular  interest  in  the  p  version,  where  curved  boundaries 
must  be  modeled  exactly  or  with  sufficient  accuracy.  Consider  a  single 
curved  element  K  in  the  mesh,  where  K  is  the  image  of  the  reference  square 
iC  =  [— 1, 1]^  under  a  smooth,  invertible  mapping  Fk-  Let  V (K)  be  the  veloc¬ 
ity  space  on  the  reference  element  (e.g.  for  method  (A),  V{K)  =  [Qp{K)]'^). 
The  usual  way  to  define  the  velocity  space  is  to  take 

V{K)  =  {u|t;  =  o  V  e  V{K)}, 

with  the  pressure  space  being  defined  similarly. 

Unfortunately,  this  definition  does  not  yield  a  combination  for  which  an 
inf-sup  condition  can  be  easily  proven,  raising  questions  about  its  stability. 
Therefore,  in  [4],  we  came  up  with  an  alternative  definition.  We  divided 
the  basis  of  V (K)  into  two  sets:  the  internal  (bubble)  basis  functions  which 
vanish  on  dK  and  the  external  basis  functions  E{K)  which  are  non-zero  on 
part  of  dK.  Then  we  defined 

I{K)  =  {v\v  =  VKiv),veI{K)} 

E{K)  =  {v\v  =  voF^\veE{k)}, 

where  Vk  represents  the  Piola  transform  from  K  onto  K.  We  took  V (K) 
to  be  the  space  spanned  by  {I{K),  E{K)},  and  defined  the  pressure  as  before. 
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With  this  definition,  the  same  stability  is  obtained  as  that  for  straight¬ 
sided  elements.  Computational  results  confirmed  that  the  method  is  optimal 
in  both  h  and  p. 

Finally,  for  triangular  elements,  we  took  the  combination  [Pp]^/Pp_2.  Ex¬ 
tensive  numerical  investigation  in  [4]  indicated  that  the  stability  behaves  like 
0{p~^)  for  this  combination.  In  [7],  we  also  theoretically  analyzed  this  ele¬ 
ment,  where  we  showed  that  it  is  stable  with  respect  to  h,  while  the  p  stability 
is  no  worse  than  0(p~^).  We  also  showed  exponential  convergence  there  for 
fiuid  fiow  problems. 

Several  numerical  experiments  are  given  in  [4].  One  of  them  concerns  a 
benchmark  problem  suggested  by  Szabo  and  Babuska,  for  the  stress  extrac¬ 
tion  around  a  circular  hole.  Figure  2  shows  this  test  problem,  and  also  shows 
a  comparison  of  the  FE  solutions  obtained  by  the  standard  p  version  FEM 
(such  as  via  STRESS  CHECK),  and  one  using  our  mixed  p  version  FEM. 
(We  are  calculating  SNS,  the  sum  of  the  normal  stresses  here.) 

The  results  in  [4]  were  also  extended  to  non-linear  problems  with  large 
displacement  and  small  strain,  when  the  Poisson  ratio  is  close  to  0.5  (e.g. 
rubber).  We  showed  that  our  elements  work  better  than  the  standard  hp 
finite  element  method.  (For  instance,  for  Poisson  ratio  0.499,  the  non-linear 
solver  in  STRESS  CHECK  fails  to  calculate  accurate  stresses,  unlike  our 
method.)  Our  future  goal  is  to  investigate  the  use  of  these  mixed  methods 
in  various  applications,  such  as  non-Newtonian  fluid  flow  (see  below)  and 
viscoelasticity. 

3.  Mixed  hp  methods  for  non-Newrtonian  flow 

We  have  applied  the  results  from  [4,5,6]  to  a  model  problem  in  CFD, 
the  upper  convected  Maxwell  fluid  (which  models  viscoelastic  flow).  This  is 
generally  solved  using  a  “three  field”  mixed  formulation  (velocity,  pressure 
and  stress).  A  question  of  crucial  importance  that  then  occurs  is  how  to 
balance  these  spaces  to  ensure  stability.  The  usual  method  that  is  used  for 
analysis  is  to  consider  the  limiting  case  of  zero  relaxation  time,  and  formulate 
spaces  that  are  stable  for  this  (linearized)  case.  Although  this  method  has 
been  used  for  h  version  spaces,  it  has  not  so  far  been  used  to  design  p  or  hp 
spaces.  In  fact,  the  computational  study  of  such  flows  by  the  p/hp  methods 
has  only  recently  been  initiated,  by  Crochet  and  Khomami. 

In  [7],  we  derived  various  choices  of  hp  spaces  for  this  problem,  and  char- 
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Angle  (Radians) 


Figure  2:  Top:  Rigid  circular  inclusion  benchmark  problem  and  mesh.  Bot¬ 
tom:  The  standard  (SFEM)  and  mixed  (MFEM)  finite  element  methods 
under  p  refinement,  u  =  .4999.  SNS  along  circular  boundary.  Degree=8. 


acterized  their  stability  and  convergence  properties.  In  particular,  we  (1) 
showed  that  the  choice  of  using  elements  for  the  (discontinuous)  stresses, 

PP  for  the  (discontinuous)  pressures  and  for  the  (continuous)  velocities 
is  completely  optimal  and  stable  in  terms  of  both  h  and  p,  and  (2)  es¬ 
tablished  exponential  hp  convergence  for  several  choices  (all  results  being  for 
the  limiting  case  of  zero  relaxation  time).  We  also  analyzed  the  hp  modified 
EVSS  method,  which  allows  us  to  take  smaller  spaces  for  the  stresses  (in 
fact,  for  the  optimal  choice  above,  we  can  now  take  only  the  reduced  space 
QP'  for  the  stresses). 

In  addition  to  formulating  these  new  combinations  of  spaces,  our  work 
also  provided  a  mathematical  validation  of  spaces  similar  to  those  used  by 
Khomami.  Here,  continuous  stresses  are  generally  used,  and  we  obtained 
estimates  for  the  following  combinations: 

For  quadrilaterals:  stresses,  P^  or  pressures  and  velocities. 

For  triangles:  PP+^  stresses,  pp~^  pressures  and  PP+^  velocities. 

Finally,  we  showed  that  for  graded  geometric  meshes,  we  obtain  exponen¬ 
tial  hp  convergence  computationally  for  a  stationary  Newtonian  flow  in  an  L- 
shaped  domain.  This  models  the  singularities  arising  also  in  non-Newtonian 
flow  at  the  reentrant  corners  in  the  so-called  4:1  contraction  problem  (see 
Figure  3). 

4.  Linearized  analysis  of  non-linear  buckling 

The  area  of  failure  prediction,  in  the  context  of  e.g.  electronic  compo¬ 
nents,  ceramics,  laminated  composites,  etc,  has  been  of  long-standing  impor¬ 
tance  to  USAF  needs.  Experimental  determination  of  loads  that  can  cause 
buckling  and  other  failures  is  both  expensive  and  unreliable.  Mathematical 
models  used  in  current  engineering  practice  are,  on  the  other  hand,  based 
on  dimensionally  reduced  descriptions  of  an  elastic  body.  The  assumptions 
necessary  for  such  dimensional  reduction  to  hold  can  often  result  in  a  large, 
unknown  modeling  error  when  non-zero  initial  stress  states  are  present,  lead¬ 
ing  to  results  such  as  unnecessary  weight  penalties  and  unexpected  failures. 

We  have  been  mathematically  validating  a  method  proposed  by  Barna 
Szabo  that  does  buckling  analysis  for  the  fully  three-dimensional  problem 
at  hand,  rather  than  some  asymptotic  (dimensionally  reduced)  limit.  This 
method  determines  buckling  loads  and  deformations  accurately  and  inex¬ 
pensively,  using  hp  methods  to  obtain  high  rates  of  convergence  even  in  the 
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presence  of  singularities. 

The  formulation  by  Szabo  reduces  to  the  question  of  finding  the  spectral 
values  of  a  non-compact  operator.  Such  problems  can  be  quite  difficult  to 
approximate  numerically,  due  to  the  presence  of  continuous  and  residual 
spectra,  and  due  to  spurious  eigenvalue  approximations.  The  first  question 
that  must  be  addressed  for  a  non-compact  operator  is  whether  the  continuous 
or  residual  spectrum  exist,  or  whether  things  are  better  behaved,  and  one  has 
only  eigenvalues  (note  that  approximate  methods  generally  are  not  succesful 
if  non-empty  continuous  or  residual  spectra  need  to  be  calculated). 

Under  this  grant,  we  established  that  in  regions  of  ellipticity  (which  are 
the  regions  of  interest),  the  spectrum  only  consists  of  real,  isolated  eigenvalues 
of  finite  multiplicity.  This  clears  the  way  for  approximation  by  numerical 
methods.  Moreover,  we  proved  that  for  every  eigenvalue  A,  there  exists  a 
sequence  of  approximate  eigenvalues  converging  to  A  when  the  finite  element 
method  is  used. 

We  are  currently  dealing  with  the  question  of  spurious  eigenvalues.  For 
certain  situations,  we  can  show  that  these  spurious  modes  will  not  pollute 
the  lowest  eigenvalue  (which  is  the  one  of  interest).  For  other  situations, 
we  have  developed  an  algorithm  that  can  detect  (and  hence  eliminate)  such 
spurious  modes.  Ongoing  research  will  concentrate  on  questions  such  as  the 
treatment  of  non-conservative  loads,  as  well  as  the  question  of  how  to  deal 
with  singular  domains.  Upon  completion  of  this  research,  we  will  have  both 
a  robust  numerical  method  suitable  for  implementation  in  the  commercial 
hp  code  STRESS  CHECK,  and  a  mathematically  rigorous  formulation  of 
buckling  of  three-dimensional  bodies,  which  incorporates  classical  engineer¬ 
ing  derivations. 

5.  p  version  space  enrichment  methods 

Space  enrichment  methods  have  traditionally  been  used  in  the  context 
of  the  h  version.  In  a  recent  paper  by  Szabo  and  Volpert,  a  p  version  space 
enrichment  method  was  shown  to  be  superior  to  the  usual  hp  version  for  con¬ 
tact  problems.  These  (non-linear)  problems  occur  in  a  variety  of  situations 
(e.g.  gears  and  machinery  of  every  kind),  and  are  related  to  wear  and  failure 
of  components. 

We  completed  the  collaboration  [8]  with  Costabel  and  Dauge  (University 
of  Rennes,  France).  In  our  paper,  we  obtained  a  theoretical  justification  of 
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the  p  version  space  enrichment  method  in  one  dimension,  which  shows  why 
the  2-d  formulation  by  Szabo  and  Volpert  works  so  well.  We  also  derived 
an  integral  equation  formulation  which  can  be  used  with  space  enrichment 
to  give  robust  exponential  convergence  when  the  singularities  lie  in  a  known 
range.  Computational  results  were  also  provided. 

6.  The  hp  approximation  of  boundziry  layers. 

Given  the  hp  capability  available  through  several  hp  codes,  we  showed 
how  the  meshes  and  degrees  should  be  designed  when  boundary  layers  are 
present.  This  formed  the  Ph.D.  thesis  of  our  student,  Christos  Xenophon- 
tos  [9],  which  built  upon  the  one-dimensional  results  obtained  in  [10].  The 
problem  Xenophontos  investigated  was 

—d^Au  +  u  =  f  in  fi,  m  =  0  on 

for  which  he  proved  exponential  convergence  in  d.  He  looked  both  at  the 
case  of  a  smooth  domain  and  the  case  of  a  square  (unsmooth)  domain  (these 
results  appearing,  respectively,  in  [11], [12]).  In  addition,the  application  of 
our  methods  to  plate  and  shell  problems  was  also  carried  out  [13],  since 
these  two-dimensional  models  have  boundary  layers  near  the  edges  as  well. 
A  key  aspect  here  was  to  choose  the  meshes  and  the  degrees  such  that  both 
boundary  layers  and  corner  singularities  are  optimally  approximated. 

Let  us  illustrate  the  main  idea  for  a  square.  In  Figure  4,  we  have  shown 
four  possible  meshes  for  a  quarter  of  a  square  plate,  of  thickness  d.  It  is 
assumed  that  polynomials  of  degree  p  are  being  used  {p  version),  while  k  is 
a  constant.  The  first  mesh  shows  the  refinement  needed  to  treat  just  the 
boundary  layer.  The  second  is  the  mesh  that  would  be  used  just  to  treat  the 
corner  singularity.  The  third  superimposes  the  first  two  meshes.  The  fourth 
is  the  mesh  that  really  should  be  used,  taking  into  account  the  nature  of 
the  singularities  at  the  corner  which  essentially  behave  like  (^)“.  Numerical 
experiments  in  [9,13]  indicate  that  whereas  the  first  mesh  is  adequate  for 
getting  a  good  error  in  the  energy  norm,  it  is  the  fourth  mesh  that  is  necessary 
if  point-wise  convergence  in  the  stresses  is  desired. 

Our  results  were  incorporated  into  the  user  guidelines  for  STRESSCHECK 
(see  the  STRESSCHECK  manual  and  the  technical  brief  by  C.  Schwab  from 
ESRD). 
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Scheme  3:  BL-corner  mesh  Scheme  4:  BL-refined  corner  mesh 


Figure  4:  Meshes  on  a  square.  Scheme  1:  Boundary  layer  mesh.  Scheme 
2:  Geometric  refinement  for  d  =  1,  with  n  =  3  layers.  Scheme  3:  Union  of 
Scheme  1  and  Scheme  2.  Scheme  4-  Union  of  Scheme  1  and  d-dependent 
geometric  refinement  with  n  =  3  layers. 
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7.  Locking-free  finite  elements  for  plates  and  shells 

Locking  is  the  phenomenon  by  which  the  numerical  approximation  of 
parameter-dependent  problems  deteriorates  for  values  of  the  parameter  close 
to  a  limiting  value. 

Joint  work  with  J.  Pitkaranta  [14]  was  motivated  by  the  locking  con¬ 
straints  found  in  shell  problems.  So  far,  satisfactory  mixed  method  tech¬ 
niques  to  handle  such  constraints  have  not  been  developed.  A  natural  strat¬ 
egy  would  be  to  extend  the  derivation  of  plate  elements  to  the  shell  set¬ 
ting.  However,  mixed  methods  for  the  Reissner-Mindlin  plate  are  derived  by 
appealing  to  a  Hemholtz  decomposition,  an  indirect  technique  without  an 
apparent  analog  for  shells. 

In  our  work,  we  developed  a  new  means  of  deriving  mixed  elements  for  the 
Reissner-Mindlin  plate,  which  does  extend  to  shells.  For  the  first  time,  the 
reduction  operators  and  polynomial  spaces  used  to  treat  the  constraints  are 
directly  derived.  Our  analysis  gives  minimal  conditions  on  plate  elements 
needed  to  prevent  locking,  and  also  to  approximate  the  boundary  layers 
present.  Using  our  analysis,  a  clear  picture  emerges  for  the  comparison  of 
currently  available  low-order  {h  version)  plate  elements,  such  as  the  MITC 
and  Arnold-Falk  ones. 

In  [15],  we  further  extended  this  work,  by  comparing  three  low-order 
plate  elements  in  their  ability  to  deal  not  only  with  locking,  but  also  bound¬ 
ary  layers.  We  found  that  the  MITC  elements  are  more  robust  in  this  regard 
than  the  so-called  Arnold-Falk  or  Arnold-Brezzi  elements.  We  also  showed 
how  these  elements  could  be  modified  at  the  boundary  to  regain  robustness. 
These  results  are  now  being  used  to  analyze  shell  problems,  where  the  inter¬ 
action  between  locking  and  boundary  layers  is  more  complicated. 

8.  Work  completed  by  C.  Schwab 

Christoph  Schwab  was  supported  under  this  grant  for  the  first  year  only, 
after  which  he  took  up  a  position  at  ETH  Zurich.  In  addition  to  collabo¬ 
rating  on  projects  3  and  6  above,  let  us  mention  two  other  projects  that  he 
completed  while  on  this  grant. 

1.  Wavelet  Galerkin  Discretization  of  Boundary  Integral  Equations  on  Sur¬ 
faces 
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In  joint  work  with  T.  von  Petersdorff  and  R.  Schneider  [16],  a  fully  dis¬ 
crete  Multiscale  (wavelet)  method  for  boundary  integral  equations  on  sur¬ 
faces  in  was  analyzed.  It  was  shown  that  sparse  stiffness  matrices  with 
0{N{logN)'^)  nonzero  entries  suffice  to  preserve  the  full  asymptotic  conver¬ 
gence  rate  of  the  scheme.  An  explicit  quadrature  strategy  was  presented, 
which  allows  the  computation  of  these  entries  with  the  necessary  accuracy 
in  0{N (log N^)  steps  (here,  N  denotes  the  number  of  degrees  of  freedom 
on  the  boundary  surface).  These  results  are  true  for  a  wide  class  of  integral 
equation  formulations  stemming  from,  e.g.,  potential  theory,  elasticity,  ex¬ 
terior  Stokes  flow,  Maxwell’s  equations  and  Helmholtz  equations  (moderate 
frequencies) . 

2.  Hierarchical  modeling  in  mechanics 

A  project  on  hierarchical  modeling  for  various  problems  in  mechanics  was 
completed,  under  partial  sponsorship  of  AFOSR.  An  important  aspect  of  this 
work  was  to  develop  estimators  that  could  be  used  to  predict  elements  where 
the  model  order  needed  to  be  increased.  The  review  article  [17]  was  written 
for  a  summer  school  in  Leicester,  England,  and  forms  one  of  six  chapters  in 
the  proceedings. 


5.  Participating  Personnel 

In  addition  to  the  Pis  Manil  Suri  and  Christoph  Schwab,  the  following  people 
were  associated  with  the  research:  ' 

Collaborators:  Barna  Szabo  (University  of  Washington,  St.  Louis,  MO), 
Ivo  Babuska  (University  of  Texas,  Austin,  TX),  Monique  Dauge  and  Martin 
Costabel  (University  of  Rennes,  Rennes,  France),  Juhani  Pitkaranta  (Helsinki 
University  of  Technology,  Helsinki,  Finland) 

Ph.D.  students:  Christos  Xenophontos,  Lt.  Col.  Lawrence  Chilton,  USAF 
and  Padmanabhan  Seshaiyar. 

Summer  salary  support  was  provided  for  P.  Seshaiyer  through  the  current 
grant. 


6.  Theses  and  journal  publications 

NOTE:  Most  publications,  including  recent  preprints,  can  be  retrieved  from 
M.  Suri’s  web  page,  http://www.math.umbc.edu/~suri. 
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1.  P.  Seshaiyer  and  M.  Suri,  “Uniform  hp  convergence  results  for  the  mor¬ 
tar  finite  element  method,”  (1997  preprint,  to  appear  in  Mathematics 
of  Computation,  1999). 

2.  P.  Seshaiyer  and  M.  Suri,  “Convergence  results  for  non-conforming  hp 
methods:  The  mortar  finite  element  method,”  Contemporary  Mathe¬ 
matics,  218,  467-473  (1998). 

3.  P.  Seshaiyer,  “Non-conforming  hp  finite  element  methods,”  Ph.D.  Dis¬ 
sertation,  University  of  Maryland  Baltimore  County,  expected  July 
1998. 

4.  L.  Chilton,  “Locking-free  mixed  hp  finte  element  methods  for  linear 
and  geometrically  nonlinear  elasticity,”  Ph.D.  Dissertation,  University 
of  Maryland  Baltimore  County,  May,  1997. 

5.  L.  Chilton  and  M.  Suri,  “On  the  selection  of  a  locking-free  hp  element 
for  elasticity  problems,” /ni.  J.  Numer.  Meth.  Eng.  40,  2045-2062 
(1997). 

6.  L.  Chilton  and  M.  Suri,  “Stable  and  optimal  mixed  hp  methods  over 
curvilinear  domains,”  (in  preparation). 

7.  C.  Schwab  and  M.  Suri,  “Mixed  hp  finite  element  methods  for  problems 
in  non-Newtonian  and  Stokes  fiows,”  (1997  preprint,  to  appear  in  a 
special  issue  of  Comp.  Meth.  Appl.  Mech.  Eng.  devoted  to  high-order 
CFD  methods,  1999). 

8.  M.  Costabel,  M.  Dauge  and  M.  Suri,  “Numerical  approximation  of  a 
singularly  perturbed  contact  problem,”  (1996  preprint,  to  appear  in 
Comp.  Meth.  Appl.  Mech.  Eng,  1998). 

9.  C.  Xenophontos,  “The  hp  version  of  the  finite  element  method  for  singu¬ 
larly  perturbed  problems,”  Ph.D.  Dissertation,  University  of  Maryland 
Baltimore  County,  April,  1996.  (thesis  can  be  retrieved  from  the  web 
address  http://www.clarkson.edu/~christos/thesis.ps) 

10.  C.  Schwab  and  M.  Suri,  “The  p  and  hp  versions  of  the  finite  element 
method  for  problems  with  boundary  layers,”  Mathematics  of  Compu¬ 
tation,  65,  1403-1429  (1996). 
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11.  C.  Xenophontos,  “The  hp  finite  element  method  for  singularly  per¬ 
turbed  problems  in  smooth  domains,”  Math  Models  and  Methods  in 
Applied  Sciences,  Vol  2,  No.  8  ,  1998. 

12.  C.  Xenophontos,  “The  hp  finite  element  method  for  singularly  per¬ 
turbed  problems  in  non-smooth  domains,”  (1997  preprint,  submitted 
to  Numerical  Methods  for  PDFs). 

13.  C.  Schwab,  M.  Suri  and  C.  Xenophontos,  “The  hp  finite  element  method 
for  problems  in  mechanics  with  boundary  layers,”  (1996  preprint,  to 
appear  in  Comp.  Meth.  Appl.  Mech.  Eng,  1998). 

14.  J.  Pitkaranta  and  M.  Suri,  “Design  principles  and  error  analysis  for 
reduced-shear  plate-bending  finite  elements,”  Numerische  Mathematik, 
75,  223-266  (1996). 

15.  J.  Pitkaranta  and  M.  Suri,  “Upper  and  lower  plate-bending  error  bounds 
for  some  plate-bending  finite  elements,”  (1998  preprint,  submitted  to 
Numerische  Mathematik,  1998). 

16.  T.  von  Petersdorff,  C.  Schwab  and  R.  Schneider,  “Multiwavelets  for 
second  kind  integral  equations,”  SIAM  J.  Numer.  Anal.,  34,  2212-2227 
(1997). 

17.  C.  Schwab,  “Hierarchic  Modeling  in  Mechanics,”  (1996  preprint,  to 
appear  as  a  survey  chapter  in  the  Proceedings  of  the  1996  Leicester 
Summer  School  in  Numerical  Mathematics,  Springer- Verlag) . 

7.  Interactions/Transitions 

Presentations  at  meetings,  conferences  and  seminars. 

Mar  ’95  Texas  A  &  M  University,  College  Station,  TX,  “The  hp  finite  ele¬ 
ment  modeling  of  plates  and  shells”  (seminar) 

Mar  ’95  University  of  Texas,  Austin,  TX,  “The  hp  finite  element  modeling 
of  plates  and  shells”  (seminar) 

Jul  ’95  ICIAM  ’95  (International  Congress  of  Industrial  and  Applied  Math¬ 
ematics),  Hamburg,  July  3-7,  1995,  “Mixed  hp  methods  for  plates  and 
shells”  (Minisymposium  lecture,  organized  minisymposium) 
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Aug  ’95  “Waves  and  Memory  in  Continua:  A  Meeting  in  Honor  of  Richard 
C.  MacCamy,”  Carnegie-Mellon  University,  Pittsburgh,  Aug  17-19,  “The 
finite  element  modeling  of  thin  structures”  (invited  plenary  lecture) 

Oct  ’95  “1995  SIAM  Annual  Meeting,”  Charlotte,  North  Carolina,  Oct  23- 
26,  “The  design  of  locking-free  reduced-shear  plate  elements”  (invited 
minisymposium  lecture) 

Dec  ’95  Courant  Institute,  New  York  University,  New  York,  NY,  “The  hp 
finite  element  modeling  of  thin  structures”  (seminar) 

Jan  ’96  Washington  University,  St.  Louis,  MO,  “The  hp  finite  element 
modeling  of  thin  domains”  (seminar) 

Mar  ’96  “Engineering  Problems:  Mathematical  Formulation,  Analysis,  and 
their  Computational  Treatment,”  University  of  Maryland  College  Park, 
Mar  21-24,  “The  hp  finite  element  modeling  of  thin  structures”  (invited 
plenary  lecture) 

Jul  ’96  “Prague  Mathematical  Conference  1996,”  Prague,  Czech  Republic, 
Jul  8-12,  “The  hp  finite  element  modeling  of  thin  structures”  (invited 
plenary  lecture) 

Jul  ’96  “Numerical  Methods  and  Computational  Mechanics  in  Science  and 
Engineering,”  Miskolc,  Hungary,  Jul  15-19,  “The  hp  approximation  of 
boundary  layers”  (invited  minisymposium  lecture) 

Aug  ’96  Ecole  Polytechnique  Federale  de  Lausanne,  Lausanne,  Switzerland, 
“The  hp  finite  element  modeling  of  thin  domains”  (seminar) 

Sep  ’96  University  of  Maryland  College  Park,  College  Park,  MD,  “The  de¬ 
sign  of  locking-free  reduced-shear  plate  elements”  (seminar) 

Oct  ’96  University  of  Texas,  Austin,  TX,  “The  design  of  locking-free  reduced- 
shear  plate  elements”  (seminar) 

Nov  ’96  University  of  Colorado,  Denver,  CO,  ‘The  hp  finite  element  mod¬ 
eling  of  thin  domains”  (seminar) 
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Apr  ’97  American  Mathematical  Society  Meeting,  College  Park,  MD,  April 
12-13,  “Uniform  hp  estimates  over  partitioned  domains”  (invited  min¬ 
isymposium  lecture) 

Apr,  ’97  Semi-annual  FEM  conference.  New  York  University,  New  York, 
NY,  “Eigenvalue  analysis  for  a  three-dimensional  buckling  model”  (con¬ 
tributed  lecture) 

Jun  ’97  “First  International  Congress  of  the  International  Society  for  Anal¬ 
ysis,  Applications  and  Computing,”  University  of  Delaware,  Newark, 
DE,  Jun  2-6,  “Mixed  hp  methods  for  elasticity  and  flow  problems” 
(invited  minisymposium  lecture) 

Aug  ’97  “US  Congress  on  Computational  Mechanics,”  San  Francisco,  CA, 
Aug  6-8,  “Eigenvalue  analysis  of  a  three-dimensional  buckling  model” 
(invited  minisymposium  lecture) 

Aug  ’97  “Tenth  International  Conference  on  Domain  Decomposition,”  Boul¬ 
der,  CO,  Aug  10-14,  “Uniform  hp  estimates  for  partitioned  domains” 
(contributed  minisymposium  lecture) 

Sep  ’97  Pennsylvania  State  University,  PA,  “Uniform  hp  estimates  over  par¬ 
titioned  domains”  (seminar) 

Oct,  97  Semi-annual  FEM  conference,  Cornell  University,  Ithaca,  NY,  Oc¬ 
tober  10-11,  1997,  “The  hp  version  for  non-Newtonian  fluid  flow  prob¬ 
lems”  (contributed  lecture) 

Transitions 

1.  Work  on  non-conforming  hp  methods  for  multistructures  validated  and 
helped  develop  an  interface  method  which  was  incorporated  in  the  commer¬ 
cial  code  MSC-NASTRAN  ,(a  product  of  the  MacNeal-Schwendler  Corpora¬ 
tion).  Our  results  were  used  in  designing  the  flnal  method  implemented. 
Contact:  Dr.  John  Schiermeier,  The  MacNeal-Schwendler  Corporation,  815 
Colorado  Blvd,  Los  Angeles,  CA  90041-1777  (phone:  213-259-3832) 

2.  Research  on  non-linear  buckling  was  initiated  in  response  to  questions 
posed  by  Prof.  Barna  Szabo  of  ESRD,  7750  Clayton  Road,  St.  Louis,  Mis¬ 
souri  63117  (phone:  314-645-1423).  This  work  is  geared  towards  formulating 
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an  implementation  on  the  commercial  code  STRESS  CHECK.  Our  work  on 
formulating  guidelines  for  hp  meshes  in  the  presence  of  boundary  layers  has 
already  been  used  in  this  code. 

8.  Inventions/patent  disclosures:  None. 

9.  Honors  and  Awards:  None 
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